Introduction

Rhythmic entrainment is the dynamic adaptation of motor action to rhythms observed in the environment. It is what musicians do, for example, when keeping a strict tempo, together. To study rhythmic entrainment in a controlled setting, music playing is often reduced to persons performing finger tappings in response to observed rhythms. Studies have analyzed such performances in terms of cross-correlation, cross-recurrence and quantification analysis, delay coupled oscillators, and other methods (see Demos & Palmer, 2023 for a general overview). However, we believe that an approach based on dynamic predictive models may shed new light on rhythmic entrainment data.

In this case study we use the drifting metronomes paradigm (Rosso et al., 2021) as a specific controlled way to study rhythmic entrainment. We show that for this specific paradigm, parameters of a dynamic predictive model can be estimated. In particular, we show that (simulated) tapping data of dyadic interactions can be predicted with the help of a set of non-linear ordinal differential equations (ODE) describing the change of phase in tapping (by subjects) and ticking (by metronomes). Estimating the parameters of these ODE is challenging due to the nonlinear nature of the interactions, sensitivity to initialization, and the complexity of the models. In this study, we show that estimation can be done with a regression model based on ODE. The estimation is implemented in the probabilistic programming language STAN.

The ODE-system and ODE-regression offers an approach for understanding entrainment at a causal level, thus forming a basis for explaining phenomena such as leader-following behavior, asymmetries (in relative phase representations), and instantaneous frequency (revealed in recurrence analysis) in terms of the physical parameters of coupling strength and phase delay.

In what follows, we don’t model the behavior of oscillators per se but we model the flow of phase in oscillators in terms of ODE, which resembles Kuramoto models for phase flow in coupled oscillator systems (e.g. Heggli et al., 2019). ODE-regression was inspired by the Bayesian modelling approaches (e.g. McElreath, 2020, Gelman et al., 2013), including some insights of how to work with STAN (and the R-package brms), see a number of blogs listed at the end of this article.

Methodological framework

The methodological framework adopted here involves three blocks, see Figure \(\ref{InfoStreamODE}\):

Overall methodology for using ODE system and ODE regression
Overall methodology for using ODE system and ODE regression

In this study, we consider only block (a) and block (b), to show that the ODE regression model is working. For block (c), see the forthcoming work by Rosso et al. (in preparation).

Objectives of this study

The primary objectives of the present study are as follows:

Implementation

To implement the model, Stan 2.26.1 is used. The calculations are performed on a dedicated server running R version 4.3.0 (2023-04-21) and R-Studio Server 2023.6.0.421 on Ubuntu 22.04.2 LTS. The hardware configuration of the server includes a dual AMD Epyc 74F3 CPU with a total of 48 cores (96 threads), 128 GB (8x16GB) ECC DDR4 RAM, and an Nvidia RTX 3090TI graphics card.

The drifting metronome paradigm

In the drifting metronomes paradigm, two human subjects engage in synchronized tapping while listening to a metronome through headphones. The two metronomes, one for each subject, are slightly detuned in tempo, so that their phase is drifting. In particular, one metronome ticked every 600 milliseconds, while the other ticked every 610 milliseconds. Initially, the metronomes start ticking together (in-phase), but over time, the ticks become wider apart, eventually reaching an anti-phase state where the ticks are closest together. Subsequently, the cycle is repeated, and the metronomes tick together again. This complete drifting cycle takes approximately 36.6 seconds. During the experiment, the subjects were instructed to follow the metronome they heard. There were two conditions: one where the subjects couldn’t see each other, and another where they could see each other.

Interestingly, Rosso et al. (in print) found that the subjects’ tapping behavior was influenced by seeing each other tapping, despite the instruction to follow their own metronome. In other words, visual coupling between the individuals led to an entrainment effect on tapping behavior. This suggests that visual cues from the partner’s tapping influenced the synchronization between the subjects, overriding the intended synchronization with the metronome they heard.

Specification of the ODE-system

The ODE-system, depicted in Figure 2, aims at capturing possible changes in phase due to coupling. It is characterized by a set of ordinary differential equations that define the instantaneous frequency in a network of compartments whose behavior will be conceived as oscillators. Accordingly, metronome 1, subject 1, metronome 2 and subject 2 are represented as oscillators \(m1,s1,m2\) and \(s2\) and their coupling is indicated by arrows that point to the coupled oscillator.

Each subject-oscillator is coupled to the own metronome-oscillator and to the other subject-oscillator. The strength of the coupling between the oscillators is denoted by the prefix \(k\). It represents the degree to which an observation, either auditory or visual, can entrain the corresponding action of tapping. The coupling strength ranges from 0 to 1, where a value of 0 indicates no coupling or influence, and a value of 1 indicates a strong coupling or influence. In other words, \(ks1m1\) indicates how strongly \(s1\) is coupled to \(m1\), or how strongly the observation of \(m1\) can affect the action of \(s1\).

Additionally, the coupling between oscillators allows for a phase delay, represented by the prefix \(d\). The phase delay specifies the time lag or phase shift between the oscillators and can range from \(-\frac{\pi}{2}\) to \(\frac{\pi}{2}\). It accounts for any temporal discrepancies or synchronization offsets between the oscillators.

Together, the ODE-system and its parameters capture the dynamics in network of coupled oscillators. Equation \(\ref{equation1}\) describe how the instantaneous frequency, or change of phase, evolves over time based on the coupling strengths and phase delays between the oscillators:

The oscillator system label{OscillatorSystem}
The oscillator system label{OscillatorSystem}

\[\begin{equation} \frac{d\theta_{i}}{dt} = \omega_{i} + \sum_{j=1}^4 K_{i,j} sin(\theta_{t,j} - \theta_{t,i} + D_{i,j}), with ~ \theta_i(0) = 0 \label{equation1} \end{equation}\]

where \(\theta\) is the phase, \(i,j\) is an index for elements in \(\left\{ m1,s1,s2,m2 \right\}\), and \(\omega\) is the instantaneous frequency, set as: \(\omega_i = \frac{2\pi*1000}{600}~[\frac{rad}{sec}]\) for \(i=1,2\) and \(\omega_i = \frac{2\pi*1000}{610}~[\frac{rad}{sec}]\) for \(i=3,4\). The matrix \(K\) contains the coupling parameters from rows (index \(i\)) to columns (index \(j\))(range: \([0,1]\)). Accordingly, \(K(2,1) = K_{s1,m1} = ks1m1\). The matrix \(D\) contains the delay parameters (range: \([-\frac{\pi}{2},\frac{\pi}{2}]\)). A phase delay between \(s1\) and \(s2\) can be defined but we decided to drop this parameter in subsequent simulations (see Discussion).

\[ K = \begin{bmatrix} &m1 &s1 &m2 &s2 \\ m1 &. &. &. &. \\ s1 &ks1m1 &. &. &1-ks1m1 \\ m2 &. &. &. &. \\ s2 &. &1-ks2m2 &ks2m2 &.\\ \end{bmatrix} \] \[ D = \begin{bmatrix} &m1 &s1 &m2 &s2 \\ m1 &. &. &. &. \\ s1 &ds1m1 &. &. &. \\ m2 &. &. &. &. \\ s2 &. &. &ds2m2 &.\\ \end{bmatrix} \]

In Equation 1, the logic for the phase flow of \(s1\) can be explained as follows (a similar logic holds for \(s2\)): The first term in the equation is the instantaneous frequency, which is inherited from \(m1\) (\(\omega_1\) = \(\omega_2\)), ensuring that \(s1\) follows the rhythm provided by \(m1\). The second term collects prediction errors weighted by coupling strengths. The first coupling strength defines the coupling of \(s1\) with \(m1\) (\(K_{s1,m1}\)), while the second coupling strength defines the coupling of \(s1\) with \(s2\) (\(K_{s1,s2}\)). However, the coupling strength \(K_{s1,s2}\) is the reverse of the coupling strength \(K_{s1,m1}\). The rationale is that \(s1\) has limited resources (between 0 and 1) for coupling so that higher coupling to one oscillator is at the cost of the coupling strength to the other oscillator, and vice verse. The prediction error \(\epsilon\) at \(t\) is then a phase difference between the observed oscillator and the own oscillator, mapped to the interval \([-1,+1]\) using a sinus function. Accordingly, the prediction errors for \(s1\) are: \(\epsilon_{t,s1m1} = sin(\theta_{t,m1} - \theta_{t,s1} + ds1m1)\) and \(\epsilon_{t,s1s2} = sin(\theta_{t,s2} - \theta_{t,s1})\). The change in \(\omega_2\) for \(s1\) at \(t\) is then determined by \(ks1m1*\epsilon_{t,s1m1}\) and \((1-ks1m1)*\epsilon_{t,s1s2}\).

Phase delays reflect the idea that \(s1\) may have a tendency to systematically tap before or after the tic of \(m1\). Recall that taps and ticks occur at times where the phase is \(2\pi\), and where the behavior of the oscillator can be compared with tap data. Anyhow, phase delays are calculated at each time instance in the equation. As they appear in the term that defines the coupling strength they will capture the difference but have otherwise little effect on the coupling strength as such.

In summary, the change in phase for \(s1\) involves its instantaneous frequency and prediction errors, influenced by coupling strengths, and phase delays.

Following the same logic, Equation \(\ref{equation1}\) describes the phase change for each oscillator. By numerical integration of the ODE, one can obtain continuous time series representing the phase dynamics of the oscillators. From these continuous time series, discrete taps can be generated by selecting the time values at multiples of \(2\pi\) radians.

ODE_regression

In ODE-regression, the objective is to estimate the ODE parameters specified in K and D, and predict the latent tapping process based on observed discrete taps. By estimating the tapping process, the ODE-regression approach aims to uncover the underlying dynamics and parameters of the tapping behavior, allowing for a better understanding of rhythmic entrainment.

The input to the ODE-regression is a sequence of phase values at a particular time, \(\theta_n = \theta_{t,i}\), with \(n=1,2,3...\) indicating a row having \(i\) as index for an oscillator (\(\left\{ m1,s1,s2,m2 \right\}\)), and \(t\) as time at which the tic/tap happens. For computational purposes, the phase is represted on a scale from \([-\pi,+\pi]\) .

The ODE_regression is then defined by a fitting based on:

\[\begin{equation} \theta_n \sim f(\eta_n, \kappa) \end{equation}\]

with \(\theta_n\) being the response at row \(n\), containing time information at phase \(2\pi\) or \(0 [rad]\). \(f\) is a von-Mises distribution, with mean \(\eta_n\), and precision \(\kappa\), realizing a mapping to the circular representation of phase in \(\theta_n\). \(\eta_n\) is the outcome of an ODE-solver \(Z\) up to the time indicated in row \(n\). Modulo \(2\pi\) (fmod) is used to prevent an undefined growing phase value. Accordingly:

\[ \eta_n = fmod(Z(t_n;\Theta),2\pi). \] The ODE-solver \(Z\) returns successive phases for each oscillator and is defined by a set of parameters \(\Theta\)

\[ \Theta = (\theta_i(0),K_{i,j}, D_{i,j}) \] with \(\theta_{m1,m2}(0) = 0\) and \(\theta_{s1,s2}(0)\) to be estimated. Obviously, parameters in K, D are also estimated. Finally, the solution of the ODE-solver up to \(t_n=T\) can be specified as \[ \\ \begin{aligned} &\theta_{t_n > 0,i} = \theta_i(0) + \int_{t>0}^{T} \left[ \omega_{i} + \sum_{j=1}^4 K_{i,j} sin(\theta_{t,j} - \theta_{t,i} + D_{i,j}) \right]~ dt &\\ &\theta_i(0) = 0&\\ \end{aligned} \]

The priors for the parameters in K and D are defined using previous available knowledge. As subjects are requested to follow their own metronome, the prior for the coupling parameter to the own metronome is high. Phase delay parameters are unknown but we assume that it is normal distributed with standard deviation 0.25. The initial phase of \(s1\) and \(s2\) has an uninformed prior. \[ \begin{aligned} & ks1m1 &\sim & normal(.9,.2) &\text{Prior for coupling s1 to m1, range: [0,1] }\\ & ks2m2 &\sim &normal(.9,.2) &\text{Prior for coupling s2 to m2, [0,1] }\\ & ds1m1 &\sim &normal(.0,.25) &\text{Prior for delay s1 to m1},[-\frac{\pi}{2},\frac{\pi}{2}]\\ & ds2m2 &\sim &normal(.0,.25) &\text{Prior for delay s2 to m2}, [-\frac{\pi}{2},\frac{\pi}{2}]\\ & \kappa &\sim &lognormal(0.1,1.0) &\text{Prior for precision} \\ & \theta_{s1}(0) &\sim &normal(.0,1.0) &\text{Prior for initial phase of s1} \\ & \theta_{s2}(0) &\sim &normal(.0,1.0) &\text{Prior for initial phase of s2} \\ \end{aligned} \]

Data generation

The data generation process consists of three steps:

These steps collectively generate a dataset from the continuous phase dynamics, as discrete metronome ticks, subject taps, and added noise, providing a realistic representation of the tapping process.

However, to simulate human tapping data it is important to constrain the data space. Constraints in coupling strengths are related to the task in the drifting metronomes paradigm, where subjects are asked to tap along with the metronome heard through their headphones. The bulk of coupling parameters may be expected having a coupling strength between about .5 and .8, but we allow for more extreme cases with low coupling to the own metronome (< .5).
Furthermore we assume for testing purposes that both subjects have the same distribution of coupling strength. According, we first generate coupling strength values for subject 1. Then we reshuffle these values and assign them to subject 2.

Constraints in phase delays can be justified by studies showing that humans tend to anticipate to the beat (see Repp and Su (2012) for a rewiew). This finding was confirmed in Rosso et al. (2023). When generating data, our goal is to sample from a distribution that roughly covers this anticipation.

Step 1: ODE phase-functions

ODE phase-functions are generated using defined parameters for coupling and phase delay. Prior knowledge informs the selection of these parameters, ensuring the generated phase-functions to align with the expected dynamics of the tapping process. For testing purposes, the coupling and phase delay parameters are specified using the distributions shown in Figure 3. This distribution allows for flexibility in capturing the range of coupling strengths and phase delays, in an attempt to generate a realistic interactions between the oscillators.

Densities for coupling parameters $ks1m1$ and $ks2m2$, and delay parameters $ds1m1$ and $ds2m2$

Densities for coupling parameters \(ks1m1\) and \(ks2m2\), and delay parameters \(ds1m1\) and \(ds2m2\)

The phase delay can be assumed to be zero, but it is often more realistic to consider a small anticipation in \(ds1m1\) and \(ds2m2\) by using a skewed normal distribution. This distribution allows for capturing the asymmetry and slight anticipation commonly observed in tapping behavior.

In Table 1, the 40 ODE-systems, generated using 4 parameters, are shown in the first 5 columns (simulation number + 4 parameters 2 \(k\), 2 \(d\)). The first 20 ODE-systems have random coupling strength and zero phase delay, the second 20 ODE-systems have the same coupling strength as the first 20 ODE-systems, and non-zero (random) phase delay:

The number of ODE-systems is purely arbitrary and conforms our practical concerns in keeping the ODE-regression calculations within a reasonable amount of time. For each ODE-system, the phase-functions are calculated at a high sampling rate of 10,000 samples per second.

Step 2: Selecting tick and tap values

For each oscillator, a snapshot of the time values from ODE phase-functions is taken at multiples \(2\pi\). The high sampling rate was to ensure an accurate extraction of discrete time values at multiples of \(2\pi\), thus marking the taps in terms of time values. The resulting data set, with rows \(n=1...N\), contains time, phases of metronomes and subjects, simulation number (= defined by the matrices) and reference oscillator, as illustrated below:

##     time    phi_m1    phi_s1    phi_m2    phi_s2 sim       names
## 1 0.6000  6.283185  6.276223  6.180182  6.187058   1 metronome 1
## 2 0.6006  6.289468  6.282494  6.186362  6.193250   1   subject 1
## 3 0.6093  6.380575  6.373420  6.275975  6.283040   1   subject 2
## 4 0.6100  6.387905  6.380736  6.283185  6.290265   1 metronome 2
## 5 1.1999 12.565323 12.542819 12.359335 12.381545   1 metronome 1
## 6 1.2021 12.588362 12.565792 12.381995 12.404271   1   subject 1

Step 3: Adding noise

Finally, some noise is added to the times selected from these ODE phase functions. The noise has the characteristic \(N(0, \sigma)\), with \(\sigma\) defined as \(0.0045\) [s]. The parameters, phase functions, and taps, of some selected ODE-systems are shown in Figure 4.

|

|

Figure 4. ODE-systems 1, 21, 6 and 26. The left column shows the relative phase, or \(s1-m1\) and \(s2-m2\). The right column shows the instantaneus periods of \(m1,s1,m2,s2\) for cycles 1, 2 and 3.

Figure 4 shows the relative phase representation in the left column and the instantaneous period representation in the right column. Consider the parameters for ODE-system 6 and 26 (see also Table 1). The coupling parameters in ODE-system 6 result in subject 2 being more strongly coupled to subject 1, causing a positive bias in the relative phase of subject 2 towards subject 1. Conversely, subject 1 is strongly coupled to its own metronome and is less influenced by subject 2, resulting in a relatively constant relative phase over time. A similar behaviour is actually reflected in the instantaneous period representation.

With the inclusion of delay lines in ODE-system 21, it can be observed that the two relative phase curves no longer overlap (compared to ODE-system 1) and that the relative phase of subject 2 is below zero, reflecting anticipation to metronome 2.

Additionally, it is important to note that the initial values of the ODE-phase function constrain it to start at \(0 ~ [rad]\). As indicated by the dashed red line, at the beginning of the second de-phasing/in-phasing cycle, the phase deviates from \(0 ~ [rad]\) but the dynamics over cycles becomes more stable (cf. cycle 2 and cycle 3). Note that the dashed red line marks the end of the first cycle.

These observations highlight the intricate dynamics and deviations from perfect synchrony in the tapping process, influenced by the coupling parameters, delay lines, and initial conditions of the ODE system.

Running STAN

The objective of the analysis is to assess the feasibility of retrieving the coupling parameters from the simulated data using ODE-regression modeling. Additionally, the fitted phase-functions, also known as posterior predictions, can be compared with the ODE phase-functions, which represent the underlying latent phase-functions of the simulated data with added noise. To facilitate this evaluation, only the second and third de-phasing/in-phasing cycles of the simulated data are considered in the dataset.

The input to stan is a list containing: number of observations (n_times), number of generated time instances (n_gen_times), initial time (time0), phase (phase), time at which the phase appears (time), generated time instances (time_gen), oscillator (name_indicator) (1=\(m1\), 2=\(s1\), 3= \(m2\), 4= \(s2\)).

## List of 7
##  $ n_times        : int 487
##  $ n_gen_times    : num 100
##  $ time0          : num 0
##  $ phase          : num [1:487] 0 0 0 0 0 0 0 0 0 0 ...
##  $ time           : num [1:487] 0 0.0127 0.0177 0.5689 0.5895 ...
##  $ time_gen       : num [1:100] 0.0001 0.3698 0.7395 1.1092 1.4789 ...
##  $ names_indicator: int [1:487] 1 4 2 2 2 1 1 4 3 3 ...

The stan function has the following specification:

  fit <- stan(
    data = StanData, # the data as input to stan
    file = filenameStancode, # the filename containing the stan code
    thin = 2,
    iter = 12000,
    warmup = 4000,
    cores = 4,
    chain = 2
  )

The call to the ODE function in STAN is:

...
vector[4] omega[n_times] = ode_rk45(do_dt, y0, time0, time, ks1m1, ks2m2, ds1m1, ds2m2) ;

...
  
 for (n in 1:n_times){
     phase[n] ~ von_mises(fmod(omega[n,names_indicator[n]] + pi(), 2*pi()) - pi(), kappa) ;
   }

The phases are retrieved with:

vector[4] phase_gen[n_gen_times] = ode_rk45(do_dt, y0, time0, time_gen, ks1m1, ks2m2, ds1m1, ds2m2) ;

Evaluating the posterior distributions

Figure 6. (a) Posterior distributions of defined and fitted coupling parameters K. (b) Error as differences between defined and fitted parameter means

An overall result of the fitted parameters is shown in 6a. The open circles show the defined coupling parameters, used for generating the data, in gray for \(ks1m1\) and in black for \(ks2m2\). The closed circles with short bar show the posterior distribution of the fitted coupling parameter, with bars representing the 95% confidence for the mean. In the ideal case, both open and cicles should overlap. In practice some do overlap, but others only partly overlap with the bar, and some do not overlap.

Figure 6b summarizes the errors. The boxplot IQR shows a median at zero and lower and upper notches at -0.06/-0.08 and +0.01/+0.01 for coupling parameters \(ks1m1\)/\(ks2m2\). It is worth mentioning that the phase delay parameters have no effect on the coupling error.

Figure 7. Overview of defined and fitted phase delay parameters D

In a similar way it is possible to show the delay parameters, see Figure 7. The figure shows that in the first 20 simulations, some fitted parameters deviate from the defined parameters. The results suggest that despite defined delays \(ds1m1\) and \(ds2m2\) are zero, solutions exist where delays have a value different from zero. In the second 20 simulations, delays are reasonable well estimated. Overall, it can be concluded that the fitting is reasonable, given random defined nature of the coupling and delay parameters.

Table 1 compares all defined and fitted parameters. The rows in bold contain a difference in coupling strength that is equal or more than .1.

Defined and fitted parameters of 40 ODE-systems. Column names have following coding convention: Defined parameters are: \(k1 = ks1m1\), \(k2 = ks2m2\), \(d1 = ds1m1\) and \(d2 = ds2m2\). Fitted parameters are: \(K1 = ks1m1\), \(K2 = ks2m2\), \(D1 = ds1m1\) and \(D2 = ds2m2\). Differences: \(Kk1 = K1 - k1\), \(Kk2 = K2 - k2\), \(Dd1 = D1 - d1\) and \(Dd2 = D2 - d1\)
defined
Fitted
Fitted - defined
sim k1 k2 d1 d2 K1 K2 D1 D2 Kk1 Kk2 Dd1 Dd2
1 0.71 0.72 0.00 0.00 0.75 0.75 0.02 0.03 0.04 0.03 0.02 0.03
2 0.76 0.61 0.00 0.00 0.73 0.61 0.03 -0.02 -0.03 0.00 0.03 -0.02
3 0.72 0.78 0.00 0.00 0.70 0.79 0.00 0.00 -0.02 0.01 0.00 0.00
4 0.88 0.46 0.00 0.00 0.88 0.46 0.03 -0.01 0.00 0.00 0.03 -0.01
5 0.80 0.71 0.00 0.00 0.86 0.73 -0.01 0.02 0.06 0.02 -0.01 0.02
6 0.67 0.29 0.00 0.00 0.77 0.32 -0.07 -0.09 0.10 0.03 -0.07 -0.09
7 0.29 0.67 0.00 0.00 0.29 0.60 -0.03 -0.07 0.00 -0.07 -0.03 -0.07
8 0.90 0.80 0.00 0.00 0.93 0.83 -0.01 0.04 0.03 0.03 -0.01 0.04
9 0.39 0.54 0.00 0.00 0.39 0.57 0.00 0.06 0.00 0.03 0.00 0.06
10 0.78 0.76 0.00 0.00 0.78 0.75 0.01 -0.02 0.00 -0.01 0.01 -0.02
11 0.71 0.82 0.00 0.00 0.69 0.84 0.01 0.02 -0.02 0.02 0.01 0.02
12 0.56 0.83 0.00 0.00 0.56 0.87 -0.01 0.02 0.00 0.04 -0.01 0.02
13 0.83 0.60 0.00 0.00 0.81 0.59 0.01 -0.02 -0.02 -0.01 0.01 -0.02
14 0.29 0.56 0.00 0.00 0.31 0.57 0.03 -0.02 0.02 0.01 0.03 -0.02
15 0.47 0.39 0.00 0.00 0.58 0.44 -0.18 0.12 0.11 0.05 -0.18 0.12
16 0.61 0.47 0.00 0.00 0.60 0.47 -0.01 -0.01 -0.01 0.00 -0.01 -0.01
17 0.82 0.29 0.00 0.00 0.97 0.40 0.04 -0.07 0.15 0.11 0.04 -0.07
18 0.54 0.71 0.00 0.00 0.53 0.68 0.04 -0.05 -0.01 -0.03 0.04 -0.05
19 0.60 0.88 0.00 0.00 0.59 0.85 0.04 0.02 -0.01 -0.03 0.04 0.02
20 0.46 0.90 0.00 0.00 0.47 0.93 -0.03 0.03 0.01 0.03 -0.03 0.03
21 0.71 0.72 -0.24 0.63 0.75 0.72 -0.25 0.62 0.04 0.00 -0.01 -0.01
22 0.76 0.61 1.01 0.35 0.78 0.59 0.98 0.34 0.02 -0.02 -0.03 -0.01
23 0.72 0.78 0.63 -0.12 0.74 0.75 0.64 -0.10 0.02 -0.03 0.01 0.02
24 0.88 0.46 0.70 1.01 0.90 0.46 0.69 0.86 0.02 0.00 -0.01 -0.15
25 0.80 0.71 0.13 0.17 0.78 0.70 0.16 0.17 -0.02 -0.01 0.03 0.00
26 0.67 0.29 0.56 0.46 0.63 0.27 0.56 0.38 -0.04 -0.02 0.00 -0.08
27 0.29 0.67 0.54 0.13 0.41 0.89 0.54 0.14 0.12 0.22 0.00 0.01
28 0.90 0.80 0.70 0.58 0.85 0.81 0.70 0.65 -0.05 0.01 0.00 0.07
29 0.39 0.54 0.18 0.51 0.38 0.54 0.20 0.54 -0.01 0.00 0.02 0.03
30 0.78 0.76 0.33 0.54 0.78 0.73 0.36 0.57 0.00 -0.03 0.03 0.03
31 0.71 0.82 0.76 0.56 0.69 0.83 0.76 0.55 -0.02 0.01 0.00 -0.01
32 0.56 0.83 0.17 0.70 0.57 0.75 0.18 0.64 0.01 -0.08 0.01 -0.06
33 0.83 0.60 0.35 0.33 0.78 0.59 0.40 0.33 -0.05 -0.01 0.05 0.00
34 0.29 0.56 0.47 0.76 0.28 0.53 0.39 0.72 -0.01 -0.03 -0.08 -0.04
35 0.47 0.39 0.29 0.18 0.56 0.44 0.12 0.26 0.09 0.05 -0.17 0.08
36 0.61 0.47 0.46 0.33 0.66 0.47 0.36 0.42 0.05 0.00 -0.10 0.09
37 0.82 0.29 -0.12 0.70 0.82 0.30 -0.11 0.50 0.00 0.01 0.01 -0.20
38 0.54 0.71 0.33 0.29 0.54 0.74 0.31 0.34 0.00 0.03 -0.02 0.05
39 0.60 0.88 0.51 0.47 0.61 0.87 0.47 0.45 0.01 -0.01 -0.04 -0.02
40 0.46 0.90 0.58 -0.24 0.47 0.89 0.51 -0.23 0.01 -0.01 -0.07 0.01

Evaluating the posterior predictions

Another evaluation is based on a more qualitative comparison of ODE phase-functions with fitted phase-functions, or posterior predictions. Below are shown some examples:

ODE 1 ODE 21 ODE 1 ODE 21 ODE 3 ODE 23 ODE 3 ODE 23

Figure 8. ODE systems 1, 21, 6 and 26. Defined phase functions (dotted lines, and circles for simulated tap data) versus fitted phase functions (full lines and gray band showing the mean and 95% confidence of the posterior predictions). The estimations are based on cycles 2 and 3.

Discussion

The error between defined and fitted coupling parameters shows a IRQ (or 50% of the ODE_systems) covering an interval of less than 2% of the coupling strength scale. This calculation is based on the mean of the posterior distributions. There is no difference in accuracy for coupling strengths with our without defined phase delay. This result suggests that the ODE_regression performans rather well for this specific ODE-system.

However, several limitations and shortcomings of the evaluation should be mentioned.

First, the definition of phase delays is purely hypothetical and even problematic when the coupling strength to the own metronome is weak. Indeed, in tapping experiments, when attention is devoted to the other subject rather than to the own metronome, phase delay to the own metronome might be less relevant and instead of anticipation it may become a real delay. Having it defined in advance is therefore somewhat exaggerated.

Second, the ODE_system assumes that in the drifting metronomes paradigm, attention is either devoted to the own metronome heard or to the other subject seen. The model of communicating vessels works rather well but one may question whether it is realistic. The communicating vessels work here between different perceptual modalities (auditory, visual) and interaction between modalities is not very well understood. Moreover, human subjects may be distracted from both the metronome and the other subject. Distraction is something that the model doesn’t capture at this moment.

Third, in the current model, the coupling strength is fixed over time. However, in realistic scenarios with human subjects, it can be assumed that the coupling strength changes over time. In the current model, this will be reflected in a large variability of the estimated posterior destribution.

Fourth, we ran several simulations which demonstrate that a phase delay between \(s1\) and \(s2\) can be estimated. However, due to the phase shift in metronomes, this phase delay will change over time and its mean will be zero when the coupling strength is 1. Typically, the parameter has considerable variance meaning that the estimation is uncertain.

The final and most important caveat is that the modelling needs to accurately reflect the physical setup. Accordingly, a slightly different paradigm for studying rhythmic entrainment (e.g. with spontaneous tapping without metronomes) will require a careful adaptation of the ODE-system used here.

Conclusion

In this study we show that parameters of a non-linear ODE-system of the Kuramoto family can be estimated using regression modelling. The posterior predictions or phase-functions generated after fitting the parameters were compared with the original phase-functions that were used to generate the discrete tapping data, and their similarity further supports the reliability of the regression model.

Overall, the study suggests that ODE-regression modelling offers valuable tool for acquiring insight about rhythmic entrainment, in particular about the parameters that allow for a causal understanding of the phase flow in dyadic interactions.

A demonstration of the ODE-regression model as a tool for analyzing human tapping data within the drifting metronome paradigm will be given in Rosso et al. (in preparation). Further research focusing on different experimental paradigms, and application of the ODE-regression model on human data has the potential to contribute to a deeper understanding of sensorimotor synchronization in music-dance (e.g. Talebi et al., 2021), motor rehabilitation (Moumdjian et al., 2022), and other forms of human expression and empowerment, including biofeedback systems (Leman, 2016).

Acknowledgment

This work is supported by my Methusalem grant at Ghent University, grant number 01M00208, BOF UGent, https://www.ugent.be/en/research/funding/bof/methusalem.

References

Shifting metronomes paradigm and sensorimotor synchronization

Demos, A. and Palmer, C. (2023). Social and nonlinear dynamics unite: musical group synchrony, Trends in Cognitive Sciences. https://doi.org/10.1016/j.tics.2023.05.005

Leman, M. (2016). The expressive moment. How interaction (with music) shapes human empowerment. Cambridge, MA: The MIT press.

Moumdjian, L, Moens, B., Van Wijmeersh, B., Leman, M., Feys, P. (2022). Application of step and beat alignment approaches and its effect on gait in progressive multiple sclerosis with severe cerebellar ataxia: A proof of concept case study. Multiple Sclerosis Journal 28 (3), 492-495.

Repp, B., Su, YH (2013). Sensorimotor synchronization: A review of recent research (2006–2012). Psychonomic Bulletin & Review 20, 403–452. https://doi.org/10.3758/s13423-012-0371-2

Rosso, M., Maes, P.J. & Leman, M. Modality-specific attractor dynamics in dyadic entrainment. Sci Rep 11, 18355 (2021). https://doi.org/10.1038/s41598-021-96054-8

Rosso et al.(in print). Embodied perspective taking enhances interpersonal synchronization. A body-swap study. iScience.

Rosso, M., Maes PJ. & Leman, M. (in preparation). (Tentative title) Estimating coupling and phase delay parameters of entrainment in a body-swap study.

Talebi, M. , Campo, A., Aarts, N. & Leman, M. (2023). Influence of musical context on sensorimotor synchronization in classical ballet solo dance. Plos one 18 (4), e0284387. https://doi.org/10.1371/journal.pone.0284387

Bayesian regression

McElreath, R. (2020). Statistical rethinking: A Bayesian course with examples in R and Stan. CRC press.

Gelman, A., Carlin, J. B., Stern, H. S., Dunson, D. B., Vehtari, A., & Rubin, D. B. (2013). Bayesian data analysis. CRC press.

Blogs

Bob Carpenter, 28 January 2018. Predator-Prey Population Dynamics: the Lotka-Volterra model in Stan. https://mc-stan.org/users/documentation/case-studies/lotka-volterra-predator-prey.html

Ben Bales & Sebastian Weber,27 July 2020. Upgrading to the new ODE interface. https://mc-stan.org/users/documentation/case-studies/convert_odes.html

Gesmann, M., and Morris, J. Hierarchical Compartmental Reserving Models. Casualty Actuarial Society, CAS Research Papers, 19 Aug. 2020, https://www.casact.org/sites/default/files/2021-02/compartmental-reserving-models-gesmannmorris0820.pdf